import numpy as np
from itertools import product
import pandas as pd

class ProductionDecisionOptimizer:
    def __init__(self):
        # 根据表2初始化数据
        self.initialize_data()
        
    def initialize_data(self):
        """根据表2初始化所有数据"""
        # 零件数据 [零件ID, 次品率, 购买单价, 检测成本]
        self.parts_data = {
            1: [0.1, 2, 1],    # [次品率, 购买单价, 检测成本]
            2: [0.1, 8, 1],
            3: [0.1, 12, 2],
            4: [0.1, 2, 1],
            5: [0.1, 8, 1],
            6: [0.1, 12, 2],
            7: [0.1, 8, 1],
            8: [0.1, 12, 2]
        }
        
        # 半成品数据 [半成品ID, 次品率, 装配成本, 检测成本, 拆解费用]
        self.semi_finished_data = {
            1: [0.1, 8, 4, 6],  # [次品率, 装配成本, 检测成本, 拆解费用]
            2: [0.1, 8, 4, 6],
            3: [0.1, 8, 4, 6]
        }
        
        # 成品数据 [次品率, 装配成本, 检测成本, 拆解费用, 市场售价, 调换损失]
        self.finished_data = [0.1, 8, 6, 10, 200, 40]
        
        # 装配关系：零件 -> 半成品
        self.assembly_relations = {
            1: 1,  # 零件1 -> 半成品1
            2: 1,  # 零件2 -> 半成品1
            3: 1,  # 零件3 -> 半成品1
            4: 2,  # 零件4 -> 半成品2
            5: 2,  # 零件5 -> 半成品2
            6: 2,  # 零件6 -> 半成品2
            7: 3,  # 零件7 -> 半成品3
            8: 3   # 零件8 -> 半成品3
        }
    
    def calculate_expected_revenue(self, x1, x2, x4):
        """计算预期收入 S - 来自检测合格而进入市场销售的成品"""
        # 根据图片中的公式：S = R ⋅ Π_{i=1 to 8}(1 - p_i x_{1i}) Π_{j=1 to 3}(1 - p_c x_{2j}) (1 - p_c x_4)
        
        # 计算零件合格的概率
        part_qual_prob = 1.0
        for i in range(8):
            pi = self.parts_data[i+1][0]  # 零件i的次品率
            # 根据公式：Π_{i=1 to 8}(1 - p_i x_{1i})
            part_qual_prob *= (1 - pi * x1[i])
        
        # 计算半成品合格的概率
        semi_qual_prob = 1.0
        for j in range(3):
            pc = self.semi_finished_data[j+1][0]  # 半成品j的次品率
            # 根据公式：Π_{j=1 to 3}(1 - p_c x_{2j})
            semi_qual_prob *= (1 - pc * x2[j])
        
        # 计算成品合格的概率
        pc_finished = self.finished_data[0]  # 成品次品率
        # 根据公式：(1 - p_c x_4)
        finished_qual_prob = 1 - pc_finished * x4
        
        # 计算最终合格产品的概率（按照图片中的公式）
        total_qual_prob = part_qual_prob * semi_qual_prob * finished_qual_prob
        
        # 预期收入
        R = self.finished_data[4]  # 市场售价
        S = R * total_qual_prob
        
        return S
    
    def calculate_expected_cost(self, x1, x2, x3, x4, x5):
        """计算预期成本 C"""
        C = 0
        
        # 1. 购买成本
        purchase_cost = sum(self.parts_data[i+1][1] for i in range(8))
        C += purchase_cost
        
        # 2. 检测成本
        # 零件检测成本
        part_inspection_cost = sum(x1[i] * self.parts_data[i+1][2] for i in range(8))
        # 半成品检测成本
        semi_inspection_cost = sum(x2[j] * self.semi_finished_data[j+1][2] for j in range(3))
        # 成品检测成本
        finished_inspection_cost = x4 * self.finished_data[2]
        
        C += part_inspection_cost + semi_inspection_cost + finished_inspection_cost
        
        # 3. 装配成本
        # 根据图片中的公式：总装配成本 = 零配件的装配成本 + 半成品的装配成本
        
        # 3.1 零配件的装配成本：Σ_{j=1 to 3} c_aj
        # 根据图片中的公式：
        # 半成品1: ca1 = (1 - p1x11)(1 - p2x12)(1 - p3x13)
        ca1 = 1.0
        for i in [0, 1, 2]:  # 零件1,2,3 -> 半成品1
            pi = self.parts_data[i+1][0]  # 零件i的次品率
            ca1 *= (1 - pi * x1[i])  # (1 - p_i x_{1i})
        
        # 半成品2: ca2 = (1 - p4x14)(1 - p5x15)(1 - p6x16)
        ca2 = 1.0
        for i in [3, 4, 5]:  # 零件4,5,6 -> 半成品2
            pi = self.parts_data[i+1][0]  # 零件i的次品率
            ca2 *= (1 - pi * x1[i])  # (1 - p_i x_{1i})
        
        # 半成品3: ca3 = (1 - p7x17)(1 - p8x18)
        ca3 = 1.0
        for i in [6, 7]:  # 零件7,8 -> 半成品3
            pi = self.parts_data[i+1][0]  # 零件i的次品率
            ca3 *= (1 - pi * x1[i])  # (1 - p_i x_{1i})
        
        # 总零配件的装配成本：Σ_{j=1 to 3} c_aj
        parts_assembly_cost = ca1 + ca2 + ca3
        
        # 3.2 半成品的装配成本：c_a ⋅ Π_{j=1 to 3}(1 - p_c x_{2j})
        ca = self.finished_data[1]  # 成品装配成本作为基础装配成本
        semi_assembly_factor = 1.0
        for j in range(3):
            pc = self.semi_finished_data[j+1][0]  # 半成品j的次品率
            semi_assembly_factor *= (1 - pc * x2[j])  # (1 - p_c x_{2j})
        
        semi_assembly_cost = ca * semi_assembly_factor
        
        # 总装配成本
        total_assembly_cost = parts_assembly_cost + semi_assembly_cost
        C += total_assembly_cost
        
        # 4. 拆解成本期望值
        # 半成品拆解成本
        semi_disassembly_cost = 0
        for j in range(3):
            pc = self.semi_finished_data[j+1][0]  # 半成品j的次品率
            cdj = self.semi_finished_data[j+1][3]  # 半成品j的拆解费用
            semi_disassembly_cost += x2[j] * pc * x3[j] * cdj
        
        # 成品拆解成本
        pc_finished = self.finished_data[0]  # 成品次品率
        c5d = self.finished_data[3]  # 成品拆解费用
        finished_disassembly_cost = x4 * pc_finished * x5 * c5d
        
        C += semi_disassembly_cost + finished_disassembly_cost
        
        # 5. 调换损失期望值
        # 调换损失是指不合格产品进入市场造成的损失
        # 根据公式：c_e ⋅ (1 - x_4) [1 - Π_{i=1 to 8}(1 - p_i x_{1i}) Π_{j=1 to 3}(1 - p_c x_{2j})]
        
        # 计算零件和半成品的合格概率
        part_qual_prob = 1.0
        for i in range(8):
            pi = self.parts_data[i+1][0]
            # 根据公式：Π_{i=1 to 8}(1 - p_i x_{1i})
            part_qual_prob *= (1 - pi * x1[i])
        
        semi_qual_prob = 1.0
        for j in range(3):
            pc = self.semi_finished_data[j+1][0]
            # 根据公式：Π_{j=1 to 3}(1 - p_c x_{2j})
            semi_qual_prob *= (1 - pc * x2[j])
        
        # 不合格产品进入市场的概率
        # 公式：1 - Π_{i=1 to 8}(1 - p_i x_{1i}) Π_{j=1 to 3}(1 - p_c x_{2j})
        defect_prob = 1 - part_qual_prob * semi_qual_prob
        
        # 调换损失：c_e ⋅ (1 - x_4) [1 - Π_{i=1 to 8}(1 - p_i x_{1i}) Π_{j=1 to 3}(1 - p_c x_{2j})]
        ce = self.finished_data[5]  # 调换损失
        replacement_loss = ce * (1 - x4) * defect_prob
        
        C += replacement_loss
        
        return C
    
    def calculate_net_profit(self, x1, x2, x3, x4, x5):
        """计算净收益 W = S - C"""
        S = self.calculate_expected_revenue(x1, x2, x4)
        C = self.calculate_expected_cost(x1, x2, x3, x4, x5)
        return S - C
    
    def optimize_brute_force(self):
        """使用暴力搜索找到最优解"""
        best_profit = float('-inf')
        best_solution = None
        best_revenue = 0
        best_cost = 0
        
        # 生成所有可能的决策组合
        # x1: 8个零件的检测决策
        # x2: 3个半成品的检测决策
        # x3: 3个半成品的拆解决策
        # x4: 成品的检测决策
        # x5: 成品的拆解决策
        
        total_combinations = 2**8 * 2**3 * 2**3 * 2 * 2
        print(f"总搜索空间: {total_combinations:,} 种组合")
        
        count = 0
        for x1_tuple in product([0, 1], repeat=8):
            for x2_tuple in product([0, 1], repeat=3):
                for x3_tuple in product([0, 1], repeat=3):
                    for x4 in [0, 1]:
                        for x5 in [0, 1]:
                            x1 = list(x1_tuple)
                            x2 = list(x2_tuple)
                            x3 = list(x3_tuple)
                            
                            profit = self.calculate_net_profit(x1, x2, x3, x4, x5)
                            
                            if profit > best_profit:
                                best_profit = profit
                                best_solution = (x1, x2, x3, x4, x5)
                                best_revenue = self.calculate_expected_revenue(x1, x2, x4)
                                best_cost = self.calculate_expected_cost(x1, x2, x3, x4, x5)
                            
                            count += 1
                            if count % 10000 == 0:
                                print(f"已搜索 {count:,} 种组合...")
        
        return best_solution, best_profit, best_revenue, best_cost
    
    def print_solution(self, solution, profit, revenue, cost):
        """打印最优解"""
        x1, x2, x3, x4, x5 = solution
        
        print("\n" + "="*60)
        print("生产决策优化结果")
        print("="*60)
        
        print(f"\n最优净收益: {profit:.2f} 元")
        print(f"预期收入: {revenue:.2f} 元")
        print(f"预期成本: {cost:.2f} 元")
        
        print(f"\n决策方案:")
        print("-" * 40)
        
        # 零件检测决策
        print("零件检测决策:")
        for i in range(8):
            status = "检测" if x1[i] == 1 else "不检测"
            cost = self.parts_data[i+1][2] if x1[i] == 1 else 0
            print(f"  零件{i+1}: {status} (成本: {cost}元)")
        
        # 半成品检测决策
        print("\n半成品检测决策:")
        for j in range(3):
            status = "检测" if x2[j] == 1 else "不检测"
            cost = self.semi_finished_data[j+1][2] if x2[j] == 1 else 0
            print(f"  半成品{j+1}: {status} (成本: {cost}元)")
        
        # 半成品拆解决策
        print("\n半成品拆解决策:")
        for j in range(3):
            status = "拆解" if x3[j] == 1 else "不拆解"
            cost = self.semi_finished_data[j+1][3] if x3[j] == 1 else 0
            print(f"  半成品{j+1}: {status} (成本: {cost}元)")
        
        # 成品检测决策
        print(f"\n成品检测决策: {'检测' if x4 == 1 else '不检测'} (成本: {self.finished_data[2] if x4 == 1 else 0}元)")
        
        # 成品拆解决策
        print(f"成品拆解决策: {'拆解' if x5 == 1 else '不拆解'} (成本: {self.finished_data[3] if x5 == 1 else 0}元)")
        
        # 成本分析
        print(f"\n成本分析:")
        print("-" * 40)
        
        # 购买成本
        purchase_cost = sum(self.parts_data[i+1][1] for i in range(8))
        print(f"购买成本: {purchase_cost}元")
        
        # 检测成本
        part_inspection = sum(x1[i] * self.parts_data[i+1][2] for i in range(8))
        semi_inspection = sum(x2[j] * self.semi_finished_data[j+1][2] for j in range(3))
        finished_inspection = x4 * self.finished_data[2]
        total_inspection = part_inspection + semi_inspection + finished_inspection
        print(f"检测成本: {total_inspection}元 (零件:{part_inspection}, 半成品:{semi_inspection}, 成品:{finished_inspection})")
        
        # 装配成本
        # 3.1 零配件的装配成本：Σ_{j=1 to 3} c_aj
        ca1 = 1.0
        for i in [0, 1, 2]:  # 零件1,2,3 -> 半成品1
            pi = self.parts_data[i+1][0]  # 零件i的次品率
            ca1 *= (1 - pi * x1[i])  # (1 - p_i x_{1i})
        
        ca2 = 1.0
        for i in [3, 4, 5]:  # 零件4,5,6 -> 半成品2
            pi = self.parts_data[i+1][0]  # 零件i的次品率
            ca2 *= (1 - pi * x1[i])  # (1 - p_i x_{1i})
        
        ca3 = 1.0
        for i in [6, 7]:  # 零件7,8 -> 半成品3
            pi = self.parts_data[i+1][0]  # 零件i的次品率
            ca3 *= (1 - pi * x1[i])  # (1 - p_i x_{1i})
        
        parts_assembly_cost = ca1 + ca2 + ca3
        
        # 3.2 半成品的装配成本：c_a ⋅ Π_{j=1 to 3}(1 - p_c x_{2j})
        ca = self.finished_data[1]  # 成品装配成本作为基础装配成本
        semi_assembly_factor = 1.0
        for j in range(3):
            pc = self.semi_finished_data[j+1][0]  # 半成品j的次品率
            semi_assembly_factor *= (1 - pc * x2[j])  # (1 - p_c x_{2j})
        
        semi_assembly_cost = ca * semi_assembly_factor
        
        total_assembly = parts_assembly_cost + semi_assembly_cost
        print(f"装配成本: {total_assembly:.2f}元 (零配件:{parts_assembly_cost:.2f}, 半成品:{semi_assembly_cost:.2f})")
        
        # 拆解成本期望值
        semi_disassembly = sum(x2[j] * self.semi_finished_data[j+1][0] * x3[j] * self.semi_finished_data[j+1][3] for j in range(3))
        finished_disassembly = x4 * self.finished_data[0] * x5 * self.finished_data[3]
        total_disassembly = semi_disassembly + finished_disassembly
        print(f"拆解成本期望值: {total_disassembly:.2f}元 (半成品:{semi_disassembly:.2f}, 成品:{finished_disassembly:.2f})")
        
        # 调换损失期望值
        # 计算零件和半成品的合格概率
        part_qual_prob = 1.0
        for i in range(8):
            pi = self.parts_data[i+1][0]
            # 根据公式：Π_{i=1 to 8}(1 - p_i x_{1i})
            part_qual_prob *= (1 - pi * x1[i])
        
        semi_qual_prob = 1.0
        for j in range(3):
            pc = self.semi_finished_data[j+1][0]
            # 根据公式：Π_{j=1 to 3}(1 - p_c x_{2j})
            semi_qual_prob *= (1 - pc * x2[j])
        
        # 不合格产品进入市场的概率
        # 公式：1 - Π_{i=1 to 8}(1 - p_i x_{1i}) Π_{j=1 to 3}(1 - p_c x_{2j})
        defect_prob = 1 - part_qual_prob * semi_qual_prob
        
        # 调换损失：c_e ⋅ (1 - x_4) [1 - Π_{i=1 to 8}(1 - p_i x_{1i}) Π_{j=1 to 3}(1 - p_c x_{2j})]
        ce = self.finished_data[5]  # 调换损失
        replacement_loss = ce * (1 - x4) * defect_prob
        
        print(f"调换损失期望值: {replacement_loss:.2f}元")

def main():
    """主函数"""
    print("生产决策优化问题求解")
    print("="*60)
    
    # 创建优化器
    optimizer = ProductionDecisionOptimizer()
    
    # 执行优化
    print("开始优化搜索...")
    solution, profit, revenue, cost = optimizer.optimize_brute_force()
    
    # 打印结果
    optimizer.print_solution(solution, profit, revenue, cost)

if __name__ == "__main__":
    main()
